# Directory and Data Import

setwd("/home/joshk/git_repositories/BrookTrout_TG/Data")
BTChallenge <- read.csv("BT Challenge_DegreeDay.csv")

{
  BTChallenge$Offspring.Acc <- as.factor(BTChallenge$Offspring.Acc)
  BTChallenge$Mother.Acc <- as.factor(BTChallenge$Mother.Acc)
  BTChallenge$Father.Acc <- as.factor(BTChallenge$Father.Acc)
  BTChallenge$Mother.ID <- as.factor(BTChallenge$Mother.ID)
  BTChallenge$Father.ID <- as.factor(BTChallenge$Father.ID)
  BTChallenge$Subject.ID <- as.factor(BTChallenge$Subject.ID)
  BTChallenge$Family <- as.factor(BTChallenge$Family)
  BTChallenge$Parent.Group <- as.factor(BTChallenge$Parent.Group)
}

options(na.action = na.omit)

## Loading in necessary packages

library("easypackages")

libraries(
  "broom", "car", "dotwhisker", "dplyr", "ggplot2", "grid",
  "gridExtra", "gtable", "itsadug", "knitr", "lme4", "lmerTest",
  "mgcv", "MuMIn", "nlme", "purrr", "rmarkdown", "splines"
)

# Adding font

library("showtext")
font_add_google(name = "Noto Sans", family = "Noto Sans")

## Additional functions

strip.lme <- function(model.name) {
  coef.names <- c(names(model.name$coefficients$fixed))
  tTab.data <- data.frame(summary(model.name)$tTable[, c(1, 2, 4)], row.names = NULL)
  tTab.data <- cbind(coef.names, tTab.data)
  return(tTab.data)
}

sig.lme <- function(model.name, sig.level = 0.1) {
  sig.data <- as.data.frame(summary(model.name)$tTable[c(which(as.numeric(paste(summary(model.name)$tTable[, 5])) <= sig.level)), ])
  return(sig.data)
}

Assigning ggplot2 themes for figures

my.theme <- theme(
  panel.grid.minor = element_blank(),
  axis.title = element_text(size = 14, family = "Noto Sans"),
  axis.text = element_text(size = 14, colour = "black", family = "Noto Sans"),
  axis.title.y = element_text(vjust = 0.5), panel.grid.major =
    element_line(colour = "grey75"), legend.title = element_text(
    size = 14,
    colour = "black", family = "Noto Sans"
  ), legend.text = element_text(
    size = 12,
    colour = "black", family = "Noto Sans"
  )
)

my.theme.dw <- theme(
  panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
  axis.title = element_text(size = 12, family = "Noto Sans"),
  axis.text = element_text(size = 12, colour = "black", family = "Noto Sans"),
  legend.title = element_text(
    size = 12,
    colour = "black", family = "Noto Sans"
  ), legend.text = element_text(
    size = 12,
    colour = "black", family = "Noto Sans"
  )
)

my.theme.diag <- theme(
  panel.grid.minor = element_blank(),
  axis.title = element_text(size = 8, family = "Noto Sans"),
  axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"),
  axis.title.y = element_text(angle = 0, vjust = 0.5), panel.grid.major =
    element_line(colour = "grey75"), legend.title = element_text(
    size = 8,
    colour = "black", family = "Noto Sans"
  ), legend.text = element_text(
    size = 8,
    colour = "black", family = "Noto Sans"
  )
)

Assessing parameter distributions by plots

Unique.Mass <- distinct(BTChallenge, Mass, .keep_all = TRUE) %>%
  select(Subject.ID, Mass, Offspring.Acc, Mother.Acc, Father.Acc, Family)

ggplot(data = Unique.Mass, aes(x = Mass)) +
  geom_histogram(colour = "black", fill = "dodgerblue2", bins = 20) +
  theme_bw() +
  my.theme.diag +
  geom_vline(
    xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
    colour = "black", linetype = "dotted"
  ) +
  xlab("Mass (g)") +
  geom_density(alpha = 0.5, aes(y = 1 / 5 * ..count..), fill = "dodgerblue2")

ggplot(data = Unique.Mass, aes(
  x = Mass, fill = as.factor(Offspring.Acc),
  colour = as.factor(Offspring.Acc)
)) +
  geom_histogram(bins = 20) +
  theme_bw() +
  my.theme.diag +
  geom_vline(
    xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
    colour = "black", linetype = "dotted"
  ) +
  xlab("Mass (g)") +
  geom_density(alpha = 0.5, aes(y = 1 / 4.2 * ..count..)) +
  scale_fill_manual(
    values =
      c("dodgerblue3", "maroon1")
  ) +
  scale_colour_manual(
    values =
      c("black", "black"), guide = FALSE
  ) +
  guides(fill = guide_legend(
    title =
      "Offspring Acclimation \n Temperature"
  ))

## Strong overlap

ggplot(data = Unique.Mass, aes(
  x = Mass, fill = as.factor(Mother.Acc),
  colour = as.factor(Mother.Acc)
)) +
  geom_histogram(bins = 20) +
  theme_bw() +
  my.theme.diag +
  geom_vline(
    xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
    colour = "black", linetype = "dotted"
  ) +
  xlab("Mass (g)") +
  geom_density(alpha = 0.5, aes(y = 2 / 7 * ..count..)) +
  scale_fill_manual(
    values =
      c("slateblue", "lightpink")
  ) +
  scale_colour_manual(
    values =
      c("black", "black"), guide = FALSE
  ) +
  guides(fill = guide_legend(
    title =
      "Maternal Acclimation \n Temperature"
  ))

ggplot(data = Unique.Mass, aes(
  x = Mass, fill = as.factor(Father.Acc),
  colour = as.factor(Father.Acc)
)) +
  geom_histogram(bins = 20) +
  theme_bw() +
  my.theme.diag +
  geom_vline(
    xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
    colour = "black", linetype = "dotted"
  ) +
  xlab("Mass (g)") +
  geom_density(alpha = 0.5, aes(y = 2 / 7 * ..count..)) +
  scale_fill_manual(
    values =
      c("seagreen3", "magenta3")
  ) +
  scale_colour_manual(
    values =
      c("black", "black"), guide = FALSE
  ) +
  guides(fill = guide_legend(
    title =
      "Paternal Acclimation \n Temperature"
  ))

## No obvious distinction between offspring from warm and cold mothers, or fathers.

Exploring metabolic responses to thermal challenge across groupings. Error bars represent +/- 95% CI.

Challenge.Sum <- BTChallenge %>%
  group_by(Challenge.Temp) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

## All offspring grouped.

ggplot(Challenge.Sum, aes(x = Challenge.Temp, y = Mean)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    color = "dodgerblue3"
  ) +
  geom_point(size = 1, colour = "maroon1") +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)")

# Robust elevation at 25 °C.

Challenge.Sum.Off <- BTChallenge %>%
  group_by(Challenge.Temp, Offspring.Acc) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

## Split by offspring acclimation temperature.

ggplot(as.data.frame(Challenge.Sum.Off), aes(
  x = Challenge.Temp,
  y = Mean, group = Offspring.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("dodgerblue", "maroon1")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature"))

# Some divergence, but does not appear systematic.

Challenge.Sum.Mat <- BTChallenge %>%
  group_by(Challenge.Temp, Mother.Acc) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

# Split by maternal acclimation temperature

ggplot(as.data.frame(Challenge.Sum.Mat), aes(
  x = Challenge.Temp,
  y = Mean, group = Mother.Acc, colour = Mother.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("slateblue", "lightpink")) +
  guides(colour = guide_legend(title = "Maternal Acclimation \n Temperature"))

# Perhaps a very subtle divergence.

Challenge.Sum.Pat <- BTChallenge %>%
  group_by(Challenge.Temp, Father.Acc) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

# Split by Paternal Acclimation Temperature.

ggplot(as.data.frame(Challenge.Sum.Pat), aes(
  x = Challenge.Temp,
  y = Mean, group = Father.Acc, colour = Father.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("seagreen3", "magenta3")) +
  guides(colour = guide_legend(title = "Paternal Acclimation \n Temperature"))

# Subtle effect at upper temperatures.

Challenge.Sum.OffMat <- BTChallenge %>%
  group_by(Challenge.Temp, Offspring.Acc, Mother.Acc) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

OffMat.Group <- c()
for (i in 1:nrow(Challenge.Sum.OffMat)) {
  if (Challenge.Sum.OffMat$Offspring.Acc[i] == 1 & Challenge.Sum.OffMat$Mother.Acc[i] == 1) {
    OffMat.Group[i] <- ("1")
  } else if (
    Challenge.Sum.OffMat$Offspring.Acc[i] == 1 & Challenge.Sum.OffMat$Mother.Acc[i] == 2) {
    OffMat.Group[i] <- ("2")
  } else if (
    Challenge.Sum.OffMat$Offspring.Acc[i] == 2 & Challenge.Sum.OffMat$Mother.Acc[i] == 1) {
    OffMat.Group[i] <- ("3")
  } else if (
    Challenge.Sum.OffMat$Offspring.Acc[i] == 2 & Challenge.Sum.OffMat$Mother.Acc[i] == 2) {
    OffMat.Group[i] <- ("4")
  }
}

Challenge.Sum.OffMat$OffMat <- OffMat.Group

# Parsing by Maternal and Offspring Acclimation Temperature

ggplot(as.data.frame(Challenge.Sum.OffMat), aes(
  x = Challenge.Temp,
  y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("black", "gold2")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

Challenge.Sum.OffPat <- BTChallenge %>%
  group_by(Challenge.Temp, Offspring.Acc, Father.Acc) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

OffPat.Group <- c()
for (i in 1:nrow(Challenge.Sum.OffPat)) {
  if (Challenge.Sum.OffPat$Offspring.Acc[i] == 1 & Challenge.Sum.OffPat$Father.Acc[i] == 1) {
    OffPat.Group[i] <- ("1")
  } else if (
    Challenge.Sum.OffPat$Offspring.Acc[i] == 1 & Challenge.Sum.OffPat$Father.Acc[i] == 2) {
    OffPat.Group[i] <- ("2")
  } else if (
    Challenge.Sum.OffPat$Offspring.Acc[i] == 2 & Challenge.Sum.OffPat$Father.Acc[i] == 1) {
    OffPat.Group[i] <- ("3")
  } else if (
    Challenge.Sum.OffPat$Offspring.Acc[i] == 2 & Challenge.Sum.OffPat$Father.Acc[i] == 2) {
    OffPat.Group[i] <- ("4")
  }
}

# Now parsing by paternal and offspring acclimation temperatures

Challenge.Sum.OffPat$OffPat <- OffPat.Group
ggplot(as.data.frame(Challenge.Sum.OffPat), aes(
  x = Challenge.Temp,
  y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

# Possible, subtle interaction.

Producing three-way “interaction” plots

Challenge.Sum.All <- BTChallenge %>%
  group_by(Challenge.Temp, Offspring.Acc, Father.Acc, Mother.Acc) %>%
  summarize(
    Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
      na.rm = TRUE
    ) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
    Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
      sqrt(length(na.omit(MO2.NMS))))
  )

OffPat.Group.All <- c()
for (i in 1:nrow(Challenge.Sum.All)) {
  if (Challenge.Sum.All$Offspring.Acc[i] == 1 &
    Challenge.Sum.All$Father.Acc[i] == 1) {
    OffPat.Group.All[i] <- ("1")
  } else if (
    Challenge.Sum.All$Offspring.Acc[i] == 1 &
      Challenge.Sum.All$Father.Acc[i] == 2) {
    OffPat.Group.All[i] <- ("2")
  } else if (
    Challenge.Sum.All$Offspring.Acc[i] == 2 &
      Challenge.Sum.All$Father.Acc[i] == 1) {
    OffPat.Group.All[i] <- ("3")
  } else if (
    Challenge.Sum.All$Offspring.Acc[i] == 2 &
      Challenge.Sum.All$Father.Acc[i] == 2) {
    OffPat.Group.All[i] <- ("4")
  }
}
Challenge.Sum.All$OffPat <- OffPat.Group.All

OffMat.Group.All <- c()
for (i in 1:nrow(Challenge.Sum.All)) {
  if (Challenge.Sum.All$Offspring.Acc[i] == 1 &
    Challenge.Sum.All$Mother.Acc[i] == 1) {
    OffMat.Group.All[i] <- ("1")
  } else if (
    Challenge.Sum.All$Offspring.Acc[i] == 1 &
      Challenge.Sum.All$Mother.Acc[i] == 2) {
    OffMat.Group.All[i] <- ("2")
  } else if (
    Challenge.Sum.All$Offspring.Acc[i] == 2 &
      Challenge.Sum.All$Mother.Acc[i] == 1) {
    OffMat.Group.All[i] <- ("3")
  } else if (
    Challenge.Sum.All$Offspring.Acc[i] == 2 &
      Challenge.Sum.All$Mother.Acc[i] == 2) {
    OffMat.Group.All[i] <- ("4")
  }
}
Challenge.Sum.All$OffMat <- OffMat.Group.All

MatPat.Group.All <- c()
for (i in 1:nrow(Challenge.Sum.All)) {
  if (Challenge.Sum.All$Mother.Acc[i] == 1 &
    Challenge.Sum.All$Father.Acc[i] == 1) {
    MatPat.Group.All[i] <- ("1")
  } else if (
    Challenge.Sum.All$Mother.Acc[i] == 1 &
      Challenge.Sum.All$Father.Acc[i] == 2) {
    MatPat.Group.All[i] <- ("2")
  } else if (
    Challenge.Sum.All$Mother.Acc[i] == 2 &
      Challenge.Sum.All$Father.Acc[i] == 1) {
    MatPat.Group.All[i] <- ("3")
  } else if (
    Challenge.Sum.All$Mother.Acc[i] == 2 &
      Challenge.Sum.All$Father.Acc[i] == 2) {
    MatPat.Group.All[i] <- ("4")
  }
}
Challenge.Sum.All$MatPat <- MatPat.Group.All

OffParError.DFacet <- ggplot(as.data.frame(Challenge.Sum.All), aes(
  x = Challenge.Temp,
  y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

## Facetting by maternal acclimation temperature first

OffParError.DFacet + facet_grid(rows = vars(Mother.Acc))

OffParError.SFacet <- ggplot(as.data.frame(Challenge.Sum.All), aes(
  x = Challenge.Temp,
  y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("slateblue", "gold2")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

## Now facetting by paternal acclimation temperature

OffParError.SFacet + facet_grid(rows = vars(Father.Acc))

ParError.OFacet <- ggplot(as.data.frame(Challenge.Sum.All), aes(
  x = Challenge.Temp,
  y = Mean, group = MatPat, shape = Mother.Acc, colour = Father.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("mediumseagreen", "mediumorchid")) +
  guides(colour = guide_legend(title = "Paternal Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

## And lastly, facetting by offspring acclimation temperature

ParError.OFacet + facet_grid(rows = vars(Offspring.Acc))

Model construction

## Building cubic regression spline with low edf (here, 4) due to low curvature. Fitting spline to global data.

plot(MO2.NMS ~ Challenge.Temp, data = na.omit(BTChallenge))

curve.model <- gam(MO2.NMS ~ s(Challenge.Temp, bs = "cs", k = 7), na.action = na.exclude, data = BTChallenge)

## Results of gam
summary(curve.model)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## MO2.NMS ~ s(Challenge.Temp, bs = "cs", k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.776467   0.005796     134   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F p-value    
## s(Challenge.Temp) 5.66      6 116.7  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.213   Deviance explained = 21.5%
## GCV = 0.08691  Scale est. = 0.086686  n = 2580
plot(curve.model)

# Appears nicely fitting.

## Pulling out predicted values

smoother <- predict(curve.model, type = "response")

With.Fit <- BTChallenge
With.Fit$MO2.NMS <- smoother
With.Fit <- With.Fit %>%
  mutate("Type" = "Fit")
Resp.Alone <- BTChallenge
Resp.Alone$Type <- "Response"
Total.Data <- rbind(With.Fit, Resp.Alone)

## Plotting the cubic regression spline

Total.Data %>%
  ggplot(aes(x = Challenge.Temp, y = MO2.NMS, fill = Type, linetype = Type)) +
  geom_point(pch = 21) +
  scale_fill_manual(values = c("black", "magenta")) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7))

## Again, appears fitting. Appending spline to dataframe

BTChallenge$bSpline <- smoother

# Constructing linear model
## Crossed random intercepts per parental ID, with crossed Offspring ID.

BTChallenge <- BTChallenge %>%
  dplyr::rename("Degree.Day" = Degree.Dat)

BTChallenge$Degree.Day <- factor(BTChallenge$Degree.Day)
Dummy <- factor(1)
Dat <- groupedData(MO2.NMS ~ 1 | Dummy, BTChallenge)

Model1 <- lme(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc, random = pdBlocked(list(pdIdent(~ 0 + Father.ID), pdIdent(~ 0 + Mother.ID), pdIdent(~ 0 + Subject.ID), pdIdent(~ 0 + Degree.Day))), na.action = na.exclude, data = Dat)

# Assuming large degrees of freedom by estimating Z values (using glmmTMB):

require("glmmTMB")

Model1TMB <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day), na.action = na.exclude, data = BTChallenge)

# Diagnosing residuals

{
  M1.Res.Alone <- ggplot(
    data = data.frame("Res" = residuals(Model1, type = "normalized")),
    aes(x = 1:nrow(Dat), y = Res)
  ) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
    xlab("Row Number") +
    ylab("Deviance Residuals")

  M1.Res.by.Fit <- ggplot(data = data.frame(
    "Res" = residuals(Model1, type = "normalized"),
    "Fit" = predict(Model1)
  ), aes(x = Fit, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
    xlab("Fitted Values") +
    ylab("Deviance Residuals")

  M1.Res.by.Resp <- ggplot(data = data.frame(
    "Res" = residuals(Model1, type = "normalized"),
    "Resp" = Dat$MO2.NMS
  ), aes(x = Resp, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
    xlab("Metabolic Rate (mg O2/hour") +
    ylab("Deviance Residuals")

  M1.QQNorm.1 <- ggplot(
    data.frame("Res" = residuals(Model1), "Fit" = predict(Model1)),
    aes(sample = Res)
  ) +
    stat_qq(colour = "gold") +
    stat_qq_line() +
    my.theme +
    theme_bw()
}

grid.arrange(M1.Res.Alone, M1.Res.by.Fit, M1.Res.by.Resp, M1.QQNorm.1, nrow = 2, top = "Metabolic Rate Residuals")

# Large df

{
  M1.Res.Alone.z <- ggplot(
    data = data.frame("Res" = residuals(Model1TMB, type = "pearson")),
    aes(x = 1:nrow(BTChallenge), y = Res)
  ) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
    xlab("Row Number") +
    ylab("Pearson Residuals")

  M1.Res.by.Fit.z <- ggplot(data = data.frame(
    "Res" = residuals(Model1TMB, type = "pearson"),
    "Fit" = predict(Model1)
  ), aes(x = Fit, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
    xlab("Fitted Values") +
    ylab("Pearson Residuals")

  M1.Res.by.Resp.z <- ggplot(data = data.frame(
    "Res" = residuals(Model1TMB, type = "pearson"),
    "Resp" = BTChallenge$MO2.NMS
  ), aes(x = Resp, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
    xlab("Metabolic Rate (mg O2/hour") +
    ylab("Pearson Residuals")

  M1.QQNorm.1.z <- ggplot(
    data.frame("Res" = residuals(Model1TMB), "Fit" = predict(Model1TMB)),
    aes(sample = Res)
  ) +
    stat_qq(colour = "gold") +
    stat_qq_line() +
    my.theme +
    theme_bw()
}

grid.arrange(M1.Res.Alone.z, M1.Res.by.Fit.z, M1.Res.by.Resp.z, M1.QQNorm.1.z, nrow = 2, top = "Metabolic Rate Residuals")

## T model: Straying from normality in upper quantiles. Fan notable in residuals by fitteds. Perhaps weight by challenge temperature. First, testing for autocorrelation.
## Z model: Bias in residuals. Very well could be a consequence of autocorrelation.

ggplot(with(
  acf(na.omit(residuals(Model1, type = "normalized")), plot = "FALSE"),
  data.frame(lag, acf)
), aes(x = lag, y = acf)) +
  geom_bar(
    stat = "identity",
    position = "identity", fill = "slateblue", colour = "black"
  ) +
  theme_bw() +
  my.theme.diag +
  ylab("Autocorrelation \n Factor")

rho <- acf(na.omit(residuals(Model1, type = "normalized")), plot = F)$acf[2]
print(rho)
## [1] 0.2207105
ggplot(with(
  acf(na.omit(residuals(Model1TMB, type = "pearson")), plot = "FALSE"),
  data.frame(lag, acf)
), aes(x = lag, y = acf)) +
  geom_bar(
    stat = "identity",
    position = "identity", fill = "slateblue", colour = "black"
  ) +
  theme_bw() +
  my.theme.diag +
  ylab("Autocorrelation \n Factor")

rho <- acf(na.omit(residuals(Model1, type = "normalized")), plot = F)$acf[2]
print(rho)
## [1] 0.2207105
# Minor autocorrelation between measures in both models model. Estimating rho.
# Faceting by paternal acclimation temperature for clearer viewing.

Dat %>%
  mutate("Res" = residuals(Model1, type = "normalized")) %>%
  ggplot(aes(x = 1:nrow(Dat), y = Res, fill = Subject.ID, colour = Subject.ID)) +
  my.theme +
  theme_bw() +
  geom_point(
    size = 1, colour = "mediumseagreen",
    alpha = 0.5, pch = 21
  ) +
  xlab("Row Number") +
  ylab("Deviance Residuals") +
  theme(legend.position = "NULL") +
  geom_line(alpha = 0.4) +
  facet_grid(
    rows =
      Dat$Father.Acc
  )

Dat %>%
  mutate("Res" = residuals(Model1TMB, type = "pearson")) %>%
  ggplot(aes(x = 1:nrow(Dat), y = Res, fill = Subject.ID, colour = Subject.ID)) +
  my.theme +
  theme_bw() +
  geom_point(
    size = 1, colour = "mediumseagreen",
    alpha = 0.5, pch = 21
  ) +
  xlab("Row Number") +
  ylab("Pearson Residuals") +
  theme(legend.position = "NULL") +
  geom_line(alpha = 0.4) +
  facet_grid(
    rows =
      Dat$Father.Acc
  )

## Similar lag to models in Penney et al (In Press). Residuals, however, appear evenly spread.

# Adjusting for autocorrelation.

Model.AR <- lme(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc, random = pdBlocked(list(pdIdent(~ 0 + Father.ID), pdIdent(~ 0 + Mother.ID), pdIdent(~ 0 + Subject.ID), pdIdent(~ 0 + Degree.Day))), correlation = corAR1(), na.action = na.exclude, data = Dat)

BTChallenge$AR_Group <- factor(with(BTChallenge, paste(Father.ID, Mother.ID, Subject.ID, sep = "_")))
BTChallenge <- BTChallenge %>%
  mutate("Time" = Challenge.Temp - 14) %>%
  mutate("Time" = factor(Time)) %>%
  arrange(AR_Group, Time) %>%
  as.data.frame(.)

ModelTMB.AR <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day) + ar1(Time + 0 | AR_Group), na.action = na.exclude, data = BTChallenge)

## Diagnostics again.
# T-model

{
  MAR.Res.Alone <- ggplot(data = data.frame("Res" = residuals(Model.AR,
    type =
      "normalized"
  )), aes(x = 1:nrow(Dat), y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
    xlab("Row Number") +
    ylab("Deviance Residuals")

  MAR.Res.by.Fit <- ggplot(data = data.frame("Res" = residuals(Model.AR,
    type =
      "normalized"
  ), "Fit" = predict(Model.AR)), aes(x = Fit, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
    xlab("Fitted Values") +
    ylab("Deviance Residuals")

  MAR.Res.by.Resp <- ggplot(data = data.frame("Res" = residuals(Model.AR,
    type =
      "normalized"
  ), "Resp" = Dat$MO2.NMS), aes(x = Resp, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
    xlab("Metabolic Rate (mg O2/hour") +
    ylab("Deviance Residuals")

  MAR.QQNorm.1 <- ggplot(
    data.frame("Res" = residuals(Model.AR), "Fit" = predict(Model.AR)),
    aes(sample = Res)
  ) +
    stat_qq(colour = "gold") +
    stat_qq_line() +
    my.theme +
    theme_bw()
}

grid.arrange(MAR.Res.Alone, MAR.Res.by.Fit, MAR.Res.by.Resp, MAR.QQNorm.1, nrow = 2, top = "Metabolic Rate Residuals - Autoregressive Correlation Structure")

# Z-model

{
  MAR.Res.Alone.z <- ggplot(data = data.frame("Res" = residuals(ModelTMB.AR,
    type =
      "response"
  )), aes(x = 1:nrow(BTChallenge), y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
    xlab("Row Number") +
    ylab("Response Residuals")

  MAR.Res.by.Fit.z <- ggplot(data = data.frame("Res" = residuals(ModelTMB.AR,
    type =
      "response"
  ), "Fit" = predict(ModelTMB.AR)), aes(x = Fit, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
    xlab("Fitted Values") +
    ylab("Response Residuals")

  MAR.Res.by.Resp.z <- ggplot(data = data.frame("Res" = residuals(ModelTMB.AR,
    type =
      "response"
  ), "Resp" = BTChallenge$MO2.NMS), aes(x = Resp, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
    xlab("Metabolic Rate (mg O2/hour") +
    ylab("Response Residuals")

  MAR.QQNorm.1.z <- ggplot(
    data.frame("Res" = residuals(ModelTMB.AR), "Fit" = predict(ModelTMB.AR)),
    aes(sample = Res)
  ) +
    stat_qq(colour = "gold") +
    stat_qq_line() +
    my.theme +
    theme_bw()
}

grid.arrange(MAR.Res.Alone.z, MAR.Res.by.Fit.z, MAR.Res.by.Resp.z, MAR.QQNorm.1.z, nrow = 2, top = "Metabolic Rate Residuals - Autoregressive Correlation Structure")

# Very similar outcome to previous models. In glmmTMB model, similarities could merely be a consequence of plotted residuals not correcting for autocorrelation. Checking

plot(residuals(ModelTMB.AR, type = "response") ~ residuals(Model1TMB, type = "response"))

plot(residuals(ModelTMB.AR, type = "pearson") ~ residuals(Model1TMB, type = "pearson"))

# Yes, this appears to be the case. Using DHARMa to assess residuals

require("DHARMa")
BTChallenge_Clean <- na.omit(BTChallenge)
BTChallenge_Clean$Time_Int <- as.integer(BTChallenge_Clean$Time)

ModelTMB.AR <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day) + ar1(Time + 0 | AR_Group), na.action = na.omit, data = BTChallenge_Clean)

dharma_resid <- simulateResiduals(fittedModel = ModelTMB.AR)
plotResiduals(dharma_resid, quantreg = T) # Yes, residuals already appear cleaner.

dharma_residAR <- recalculateResiduals(dharma_resid, group = BTChallenge_Clean$Time)
testTemporalAutocorrelation(simulationOutput = dharma_residAR, time = sort(unique(BTChallenge_Clean$Time)))

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.7116, p-value = 0.5655
## alternative hypothesis: true autocorrelation is not 0
ggplot(with(
  acf(na.omit(residuals(Model.AR, type = "normalized")), plot = "FALSE"),
  data.frame(lag, acf)
), aes(x = lag, y = acf)) +
  geom_bar(
    stat = "identity",
    position = "identity", fill = "slateblue", colour = "black"
  ) +
  theme_bw() +
  my.theme.diag +
  ylab("Autocorrelation \n Factor")

## Yes, autocorrelation and residuals do appear okay in both models.
# Again, faceting by paternal acclimation temperature for clearer viewing.

Dat %>%
  mutate("Res" = residuals(Model.AR, type = "normalized")) %>%
  ggplot(aes(x = 1:nrow(Dat), y = Res, fill = Subject.ID, colour = Subject.ID)) +
  my.theme +
  theme_bw() +
  geom_point(
    size = 1, colour = "mediumseagreen",
    alpha = 0.5, pch = 21
  ) +
  xlab("Row Number") +
  ylab("Deviance Residuals") +
  theme(legend.position = "NULL") +
  geom_line(alpha = 0.4) +
  facet_grid(
    rows =
      Dat$Father.Acc
  )

BTChallenge_Clean %>%
  mutate("Res" = dharma_resid$scaledResiduals) %>%
  ggplot(aes(x = 1:nrow(BTChallenge_Clean), y = Res, fill = Subject.ID, colour = Subject.ID)) +
  my.theme +
  theme_bw() +
  geom_point(
    size = 1, colour = "mediumseagreen",
    alpha = 0.5, pch = 21
  ) +
  xlab("Row Number") +
  ylab("Deviance Residuals") +
  theme(legend.position = "NULL") +
  geom_line(alpha = 0.4) +
  facet_grid(
    rows =
      BTChallenge_Clean$Father.Acc
  )

# Good.

## Addressing tail in residuals

Dat %>%
  mutate("Res" = residuals(Model1, type = "normalized"), "Fit" = predict(Model.AR)) %>%
  ggplot(aes(sample = Res, fill = Challenge.Temp)) +
  geom_qq(pch = 21, size = 3, alpha = 0.5, aes(fill = factor(Dat$Challenge.Temp))) +
  stat_qq_line() +
  my.theme +
  theme_bw()

G_Col <- c(Cold = "dodgerblue3", Warm = "gold2", Hot = "firebrick")
Dat %>%
  mutate("Res" = abs(residuals(Model1, type = "normalized"))) %>%
  mutate("Temp_Group" = ifelse(Challenge.Temp < 20, "Cold",
    ifelse(Challenge.Temp > 25, "Hot", "Warm")
  )) %>%
  mutate("Temp_Group" = factor(Temp_Group, levels = c("Cold", "Warm", "Hot"))) %>%
  ggplot(aes(x = Res, fill = Temp_Group)) +
  geom_density(alpha = 0.4, colour = "black", adjust = 2) +
  scale_fill_manual(values = G_Col) +
  theme_classic()

BTChallenge_Clean %>%
  mutate("Res" = dharma_resid$scaledResiduals) %>%
  mutate("Temp_Group" = ifelse(Challenge.Temp < 20, "Cold",
    ifelse(Challenge.Temp > 25, "Hot", "Warm")
  )) %>%
  mutate("Temp_Group" = factor(Temp_Group, levels = c("Cold", "Warm", "Hot"))) %>%
  ggplot(aes(x = Res, fill = Temp_Group)) +
  geom_density(alpha = 0.4, colour = "black", adjust = 2) +
  scale_fill_manual(values = G_Col) +
  theme_classic()

# Slightly disproportionate pulling in cold and warm groups. Weighting by challenge temperature

Model.AR2 <- lme(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc, random = pdBlocked(list(pdIdent(~ 0 + Father.ID), pdIdent(~ 0 + Mother.ID), pdIdent(~ 0 + Subject.ID), pdIdent(~ 0 + Degree.Day))), correlation = corAR1(), weights = ~ 1 / Challenge.Temp, na.action = na.exclude, data = Dat)

ModelTMB2.AR <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day) + ar1(Time + 0 | AR_Group), weights = 1 / Challenge.Temp, na.action = na.omit, data = BTChallenge_Clean)

# glmmTMB model failing to converge. Assessing residuals of lme model.

{
  MAR2.Res.Alone <- ggplot(data = data.frame("Res" = residuals(Model.AR2,
    type =
      "normalized"
  )), aes(x = 1:nrow(Dat), y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
    xlab("Row Number") +
    ylab("Deviance Residuals")

  MAR2.Res.by.Fit <- ggplot(data = data.frame("Res" = residuals(Model.AR2,
    type =
      "normalized"
  ), "Fit" = predict(Model.AR2)), aes(x = Fit, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
    xlab("Fitted Values") +
    ylab("Deviance Residuals")

  MAR2.Res.by.Resp <- ggplot(data = data.frame("Res" = residuals(Model.AR2,
    type =
      "normalized"
  ), "Resp" = Dat$MO2.NMS), aes(x = Resp, y = Res)) +
    my.theme +
    theme_bw() +
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
    xlab("Metabolic Rate (mg O2/hour") +
    ylab("Deviance Residuals")

  MAR2.QQNorm <- ggplot(data.frame("Res" = residuals(Model.AR2, type = "normalized"), "Fit" = predict(Model.AR2)), aes(sample = Res)) +
    stat_qq(colour = "gold") +
    stat_qq_line() +
    my.theme +
    theme_bw()
}

grid.arrange(MAR2.Res.Alone, MAR2.Res.by.Fit, MAR2.Res.by.Resp, MAR2.QQNorm, nrow = 2, top = "Metabolic Rate Residuals - Autoregressive Correlation Structure and Weighted by Challenge Temp.")

## Tail persisting, but subtly. Plotting histogram of residuals.

Dat %>%
  mutate("Res" = residuals(Model.AR2, type = "normalized")) %>%
  ggplot(aes(x = Res)) +
  geom_histogram(fill = "magenta", colour = "black", size = 1, bins = 20) +
  xlab("Normalized Residuals")

## Appears quite nice, though a few notable outliers. Identifying outliers.

Dat <- Dat %>%
  mutate("Res" = residuals(Model.AR2, type = "normalized"))

Grouped <- rbind(subset(Dat, Res < with(na.omit(Dat), mean(Res) - 3 * sd(Res))) %>%
  mutate("Res.Group" = "Low.Out"), subset(Dat, Res > with(na.omit(Dat), mean(Res) + 3 * sd(Res))) %>%
  mutate("Res.Group" = "High.Out"), subset(Dat, Res >= with(na.omit(Dat), mean(Res) - 3 * sd(Res)) & Res <= with(na.omit(Dat), mean(Res) + 3 * sd(Res))) %>%
  mutate("Res.Group" = "Normal"))

G1 <- data.frame(subset(Grouped, Res.Group == "High.Out") %>%
  count(Offspring.Acc), "Group" = "Offspring.Acc")
colnames(G1) <- c("Group.Number", "Count", "Group")

G2 <- data.frame(subset(Grouped, Res.Group == "High.Out") %>%
  count(Father.Acc), "Group" = "Father.Acc")
colnames(G2) <- c("Group.Number", "Count", "Group")

G3 <- data.frame(subset(Grouped, Res.Group == "High.Out") %>%
  count(Mother.Acc), "Group" = "Mother.Acc")
colnames(G3) <- c("Group.Number", "Count", "Group")

rbind(G1, G2, G3)
##   Group.Number Count         Group
## 1            1     5 Offspring.Acc
## 2            2    10 Offspring.Acc
## 3            1    11    Father.Acc
## 4            2     4    Father.Acc
## 5            1     5    Mother.Acc
## 6            2    10    Mother.Acc
## Possible interactive assortment. Checking below.

subset(Grouped, Res.Group == "High.Out") %>%
  count(Father.Acc:Mother.Acc)
## Grouped Data: MO2.NMS ~ 1 | Dummy
##   Father.Acc:Mother.Acc n
## 1                   1:1 3
## 2                   1:2 8
## 3                   2:1 2
## 4                   2:2 2
## Not distinct. Checking trends accross acclimation temperature

subset(Grouped, Res.Group == "High.Out") %>%
  count(Challenge.Temp)
## Grouped Data: MO2.NMS ~ 1 | Dummy
##   Challenge.Temp n
## 1             18 1
## 2             19 2
## 3             20 1
## 4             21 1
## 5             23 1
## 6             25 1
## 7             26 3
## 8             27 3
## 9             28 2
## Not clear either. Residual trends do not appear systematic. Continuing with previous model as final.

Plotting Predicted Values

Dat$Pred <- predict(Model.AR)

Challenge.Sum.Pred <- Dat %>%
  group_by(Challenge.Temp, Offspring.Acc, Father.Acc, Mother.Acc) %>%
  summarize(
    Mean = mean(Pred, na.rm = TRUE), Upper.CI = mean(Pred,
      na.rm = TRUE
    ) + 1.96 * (sd(Pred, na.rm = TRUE) / sqrt(length(na.omit(Pred)))),
    Lower.CI = mean(Pred, na.rm = TRUE) - 1.96 * (sd(Pred, na.rm = TRUE) /
      sqrt(length(na.omit(Pred))))
  )

OffPat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
  if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
    Challenge.Sum.Pred$Father.Acc[i] == 1) {
    OffPat.Group.Pred[i] <- ("1")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
      Challenge.Sum.Pred$Father.Acc[i] == 2) {
    OffPat.Group.Pred[i] <- ("2")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Father.Acc[i] == 1) {
    OffPat.Group.Pred[i] <- ("3")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Father.Acc[i] == 2) {
    OffPat.Group.Pred[i] <- ("4")
  }
}
Challenge.Sum.Pred$OffPat <- OffPat.Group.Pred

# Faceting by maternal acclimation temperature.

PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
  x = Challenge.Temp,
  y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

PredError + facet_grid(rows = vars(Mother.Acc))

OffMat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
  if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
    Challenge.Sum.Pred$Mother.Acc[i] == 1) {
    OffMat.Group.Pred[i] <- ("1")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
      Challenge.Sum.Pred$Mother.Acc[i] == 2) {
    OffMat.Group.Pred[i] <- ("2")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Mother.Acc[i] == 1) {
    OffMat.Group.Pred[i] <- ("3")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Mother.Acc[i] == 2) {
    OffMat.Group.Pred[i] <- ("4")
  }
}

Challenge.Sum.Pred$OffMat <- OffMat.Group.Pred

# And faceting by paternal acclimation temperature

PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
  x = Challenge.Temp,
  y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("cornflowerblue", "gold2")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

PredError + facet_grid(rows = vars(Father.Acc))

## Striking model fit. Checking glmmTMB model

BTChallenge_Clean$Pred <- predict(ModelTMB.AR)

Challenge.Sum.Pred <- BTChallenge_Clean %>%
  group_by(Challenge.Temp, Offspring.Acc, Father.Acc, Mother.Acc) %>%
  summarize(
    Mean = mean(Pred, na.rm = TRUE), Upper.CI = mean(Pred,
      na.rm = TRUE
    ) + 1.96 * (sd(Pred, na.rm = TRUE) / sqrt(length(na.omit(Pred)))),
    Lower.CI = mean(Pred, na.rm = TRUE) - 1.96 * (sd(Pred, na.rm = TRUE) /
      sqrt(length(na.omit(Pred))))
  )

OffPat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
  if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
    Challenge.Sum.Pred$Father.Acc[i] == 1) {
    OffPat.Group.Pred[i] <- ("1")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
      Challenge.Sum.Pred$Father.Acc[i] == 2) {
    OffPat.Group.Pred[i] <- ("2")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Father.Acc[i] == 1) {
    OffPat.Group.Pred[i] <- ("3")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Father.Acc[i] == 2) {
    OffPat.Group.Pred[i] <- ("4")
  }
}
Challenge.Sum.Pred$OffPat <- OffPat.Group.Pred

# Faceting by maternal acclimation temperature again.

PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
  x = Challenge.Temp,
  y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

PredError + facet_grid(rows = vars(Mother.Acc))

OffMat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
  if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
    Challenge.Sum.Pred$Mother.Acc[i] == 1) {
    OffMat.Group.Pred[i] <- ("1")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
      Challenge.Sum.Pred$Mother.Acc[i] == 2) {
    OffMat.Group.Pred[i] <- ("2")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Mother.Acc[i] == 1) {
    OffMat.Group.Pred[i] <- ("3")
  } else if (
    Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
      Challenge.Sum.Pred$Mother.Acc[i] == 2) {
    OffMat.Group.Pred[i] <- ("4")
  }
}

Challenge.Sum.Pred$OffMat <- OffMat.Group.Pred

# Paternal acclimation temperature?

PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
  x = Challenge.Temp,
  y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
  theme_bw() +
  my.theme +
  geom_errorbar(
    mapping = aes(
      x = Challenge.Temp,
      ymin = c(Lower.CI), ymax = c(Upper.CI)
    ), width = 0.2, size = 1,
    position = position_dodge(width = 0.5)
  ) +
  geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
  xlab("Challenge Temperature \n (°C)") +
  ylab("Metabolic Rate (mg O2/h)") +
  theme(legend.position = "top") +
  scale_colour_manual(values = c("cornflowerblue", "gold2")) +
  guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
  guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
  theme(legend.title = element_text(size = 10))

PredError + facet_grid(rows = vars(Father.Acc))

# Pulling out degrees of freedom using Satterthwaite method (in lmerTest)

DF.Model <- lmer(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day), na.action = na.exclude, data = Dat)

DF <- data.frame(summary(DF.Model)$coefficients)$df
t_vals <- data.frame(summary(Model.AR)$tTab)$t.value

p_vals <- c()
for (i in 1:length(t_vals)) {
  if (t_vals[i] > 0) {
    p_vals[i] <- 2 * (pt(t_vals[i], df = DF[i], lower = FALSE))
  } else if (t_vals[i] < 0) {
    p_vals[i] <- 2 * (pt(t_vals[i], df = DF[i], lower = TRUE))
  }
}

coef_output <- data.frame(summary(Model.AR)$tTab)
coef_output$p.value <- p_vals
coef_output$DF <- DF

# Printing full model output

summary(Model.AR)
## Linear mixed-effects model fit by REML
##   Data: Dat 
##        AIC       BIC   logLik
##   -797.621 -663.1764 421.8105
## 
## Random effects:
##  Composite Structure: Blocked
## 
##  Block 1: Father.ID1, Father.ID2, Father.ID3, Father.ID4
##  Formula: ~0 + Father.ID | Dummy
##  Structure: Multiple of an Identity
##         Father.ID1 Father.ID2 Father.ID3 Father.ID4
## StdDev: 0.07769035 0.07769035 0.07769035 0.07769035
## 
##  Block 2: Mother.ID1, Mother.ID2, Mother.ID3, Mother.ID4
##  Formula: ~0 + Mother.ID | Dummy
##  Structure: Multiple of an Identity
##          Mother.ID1  Mother.ID2  Mother.ID3  Mother.ID4
## StdDev: 0.005968334 0.005968334 0.005968334 0.005968334
## 
##  Block 3: Subject.ID1, Subject.ID2, Subject.ID3, Subject.ID4, Subject.ID5, Subject.ID6, Subject.ID7, Subject.ID8, Subject.ID9, Subject.ID10, Subject.ID11, Subject.ID12, Subject.ID13, Subject.ID14, Subject.ID15, Subject.ID16, Subject.ID17, Subject.ID18, Subject.ID19, Subject.ID20, Subject.ID21, Subject.ID22, Subject.ID23, Subject.ID24, Subject.ID25, Subject.ID26, Subject.ID27, Subject.ID28, Subject.ID29, Subject.ID30, Subject.ID31, Subject.ID32, Subject.ID33, Subject.ID34, Subject.ID35, Subject.ID36, Subject.ID37, Subject.ID38, Subject.ID39, Subject.ID40, Subject.ID41, Subject.ID42, Subject.ID43, Subject.ID44, Subject.ID45, Subject.ID46, Subject.ID47, Subject.ID48, Subject.ID49, Subject.ID50, Subject.ID51, Subject.ID52, Subject.ID53, Subject.ID54, Subject.ID55, Subject.ID56, Subject.ID57, Subject.ID58, Subject.ID59, Subject.ID60, Subject.ID61, Subject.ID62, Subject.ID63, Subject.ID64, Subject.ID65, Subject.ID66, Subject.ID67, Subject.ID68, Subject.ID69, Subject.ID70, Subject.ID71, Subject.ID72, Subject.ID73, Subject.ID74, Subject.ID75, Subject.ID76, Subject.ID77, Subject.ID78, Subject.ID79, Subject.ID80, Subject.ID81, Subject.ID82, Subject.ID83, Subject.ID84, Subject.ID85, Subject.ID86, Subject.ID87, Subject.ID88, Subject.ID89, Subject.ID90, Subject.ID91, Subject.ID92, Subject.ID93, Subject.ID94, Subject.ID95, Subject.ID96, Subject.ID97, Subject.ID98, Subject.ID99, Subject.ID100, Subject.ID101, Subject.ID102, Subject.ID103, Subject.ID104, Subject.ID105, Subject.ID106, Subject.ID107, Subject.ID108, Subject.ID109, Subject.ID110, Subject.ID111, Subject.ID112, Subject.ID113, Subject.ID114, Subject.ID115, Subject.ID116, Subject.ID117, Subject.ID118, Subject.ID119, Subject.ID120, Subject.ID121, Subject.ID122, Subject.ID123, Subject.ID124, Subject.ID125, Subject.ID126, Subject.ID127, Subject.ID128, Subject.ID129, Subject.ID130, Subject.ID131, Subject.ID132, Subject.ID133, Subject.ID134, Subject.ID135, Subject.ID136, Subject.ID137, Subject.ID138, Subject.ID139, Subject.ID140, Subject.ID141, Subject.ID142, Subject.ID143, Subject.ID144, Subject.ID145, Subject.ID146, Subject.ID147, Subject.ID148, Subject.ID149, Subject.ID150, Subject.ID151, Subject.ID152, Subject.ID153, Subject.ID154, Subject.ID155, Subject.ID156, Subject.ID157, Subject.ID158, Subject.ID159, Subject.ID160, Subject.ID161, Subject.ID162, Subject.ID163, Subject.ID164, Subject.ID165, Subject.ID166, Subject.ID167, Subject.ID168, Subject.ID169, Subject.ID170, Subject.ID171, Subject.ID172, Subject.ID173, Subject.ID174, Subject.ID175, Subject.ID176, Subject.ID177, Subject.ID178, Subject.ID179, Subject.ID180, Subject.ID181, Subject.ID182, Subject.ID183, Subject.ID184, Subject.ID185, Subject.ID186, Subject.ID187, Subject.ID188, Subject.ID189, Subject.ID190, Subject.ID191, Subject.ID192, Subject.ID193, Subject.ID194, Subject.ID195, Subject.ID196, Subject.ID197, Subject.ID198, Subject.ID199, Subject.ID200, Subject.ID201, Subject.ID202, Subject.ID203, Subject.ID204, Subject.ID205, Subject.ID206, Subject.ID207, Subject.ID208, Subject.ID209, Subject.ID210, Subject.ID211, Subject.ID212, Subject.ID213, Subject.ID214, Subject.ID215, Subject.ID216, Subject.ID217, Subject.ID218, Subject.ID219, Subject.ID220, Subject.ID221, Subject.ID222, Subject.ID223, Subject.ID224, Subject.ID225, Subject.ID227, Subject.ID228, Subject.ID229, Subject.ID230
##  Formula: ~0 + Subject.ID | Dummy
##  Structure: Multiple of an Identity
##         Subject.ID1 Subject.ID2 Subject.ID3 Subject.ID4 Subject.ID5 Subject.ID6
## StdDev:   0.1471092   0.1471092   0.1471092   0.1471092   0.1471092   0.1471092
##         Subject.ID7 Subject.ID8 Subject.ID9 Subject.ID10 Subject.ID11
## StdDev:   0.1471092   0.1471092   0.1471092    0.1471092    0.1471092
##         Subject.ID12 Subject.ID13 Subject.ID14 Subject.ID15 Subject.ID16
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID17 Subject.ID18 Subject.ID19 Subject.ID20 Subject.ID21
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID22 Subject.ID23 Subject.ID24 Subject.ID25 Subject.ID26
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID27 Subject.ID28 Subject.ID29 Subject.ID30 Subject.ID31
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID32 Subject.ID33 Subject.ID34 Subject.ID35 Subject.ID36
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID37 Subject.ID38 Subject.ID39 Subject.ID40 Subject.ID41
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID42 Subject.ID43 Subject.ID44 Subject.ID45 Subject.ID46
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID47 Subject.ID48 Subject.ID49 Subject.ID50 Subject.ID51
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID52 Subject.ID53 Subject.ID54 Subject.ID55 Subject.ID56
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID57 Subject.ID58 Subject.ID59 Subject.ID60 Subject.ID61
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID62 Subject.ID63 Subject.ID64 Subject.ID65 Subject.ID66
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID67 Subject.ID68 Subject.ID69 Subject.ID70 Subject.ID71
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID72 Subject.ID73 Subject.ID74 Subject.ID75 Subject.ID76
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID77 Subject.ID78 Subject.ID79 Subject.ID80 Subject.ID81
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID82 Subject.ID83 Subject.ID84 Subject.ID85 Subject.ID86
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID87 Subject.ID88 Subject.ID89 Subject.ID90 Subject.ID91
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID92 Subject.ID93 Subject.ID94 Subject.ID95 Subject.ID96
## StdDev:    0.1471092    0.1471092    0.1471092    0.1471092    0.1471092
##         Subject.ID97 Subject.ID98 Subject.ID99 Subject.ID100 Subject.ID101
## StdDev:    0.1471092    0.1471092    0.1471092     0.1471092     0.1471092
##         Subject.ID102 Subject.ID103 Subject.ID104 Subject.ID105 Subject.ID106
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID107 Subject.ID108 Subject.ID109 Subject.ID110 Subject.ID111
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID112 Subject.ID113 Subject.ID114 Subject.ID115 Subject.ID116
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID117 Subject.ID118 Subject.ID119 Subject.ID120 Subject.ID121
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID122 Subject.ID123 Subject.ID124 Subject.ID125 Subject.ID126
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID127 Subject.ID128 Subject.ID129 Subject.ID130 Subject.ID131
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID132 Subject.ID133 Subject.ID134 Subject.ID135 Subject.ID136
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID137 Subject.ID138 Subject.ID139 Subject.ID140 Subject.ID141
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID142 Subject.ID143 Subject.ID144 Subject.ID145 Subject.ID146
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID147 Subject.ID148 Subject.ID149 Subject.ID150 Subject.ID151
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID152 Subject.ID153 Subject.ID154 Subject.ID155 Subject.ID156
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID157 Subject.ID158 Subject.ID159 Subject.ID160 Subject.ID161
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID162 Subject.ID163 Subject.ID164 Subject.ID165 Subject.ID166
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID167 Subject.ID168 Subject.ID169 Subject.ID170 Subject.ID171
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID172 Subject.ID173 Subject.ID174 Subject.ID175 Subject.ID176
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID177 Subject.ID178 Subject.ID179 Subject.ID180 Subject.ID181
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID182 Subject.ID183 Subject.ID184 Subject.ID185 Subject.ID186
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID187 Subject.ID188 Subject.ID189 Subject.ID190 Subject.ID191
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID192 Subject.ID193 Subject.ID194 Subject.ID195 Subject.ID196
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID197 Subject.ID198 Subject.ID199 Subject.ID200 Subject.ID201
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID202 Subject.ID203 Subject.ID204 Subject.ID205 Subject.ID206
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID207 Subject.ID208 Subject.ID209 Subject.ID210 Subject.ID211
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID212 Subject.ID213 Subject.ID214 Subject.ID215 Subject.ID216
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID217 Subject.ID218 Subject.ID219 Subject.ID220 Subject.ID221
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID222 Subject.ID223 Subject.ID224 Subject.ID225 Subject.ID227
## StdDev:     0.1471092     0.1471092     0.1471092     0.1471092     0.1471092
##         Subject.ID228 Subject.ID229 Subject.ID230
## StdDev:     0.1471092     0.1471092     0.1471092
## 
##  Block 4: Degree.Day15, Degree.Day30, Degree.Day45, Degree.Day60, Degree.Day75, Degree.Day90, Degree.Day105, Degree.Day120, Degree.Day135, Degree.Day150, Degree.Day165, Degree.Day180, Degree.Day195, Degree.Day210, Degree.Day285, Degree.Day304, Degree.Day323, Degree.Day342, Degree.Day361, Degree.Day380, Degree.Day399, Degree.Day418, Degree.Day437, Degree.Day456, Degree.Day475, Degree.Day494, Degree.Day513, Degree.Day532, Degree.Day551, Degree.Day570
##  Formula: ~0 + Degree.Day | Dummy
##  Structure: Multiple of an Identity
##         Degree.Day15 Degree.Day30 Degree.Day45 Degree.Day60 Degree.Day75
## StdDev:   0.06349814   0.06349814   0.06349814   0.06349814   0.06349814
##         Degree.Day90 Degree.Day105 Degree.Day120 Degree.Day135 Degree.Day150
## StdDev:   0.06349814    0.06349814    0.06349814    0.06349814    0.06349814
##         Degree.Day165 Degree.Day180 Degree.Day195 Degree.Day210 Degree.Day285
## StdDev:    0.06349814    0.06349814    0.06349814    0.06349814    0.06349814
##         Degree.Day304 Degree.Day323 Degree.Day342 Degree.Day361 Degree.Day380
## StdDev:    0.06349814    0.06349814    0.06349814    0.06349814    0.06349814
##         Degree.Day399 Degree.Day418 Degree.Day437 Degree.Day456 Degree.Day475
## StdDev:    0.06349814    0.06349814    0.06349814    0.06349814    0.06349814
##         Degree.Day494 Degree.Day513 Degree.Day532 Degree.Day551 Degree.Day570
## StdDev:    0.06349814    0.06349814    0.06349814    0.06349814    0.06349814
##          Residual
## StdDev: 0.2026764
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | Dummy 
##  Parameter estimate(s):
##       Phi 
## 0.3510209 
## Fixed effects:  MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc 
##                                                     Value  Std.Error   DF
## (Intercept)                                    -0.8447864 0.11289296 2554
## Mass                                            0.2049529 0.01670045 2554
## bSpline                                         1.1932080 0.09159720 2554
## Offspring.Acc2                                  0.0173138 0.11502032 2554
## Father.Acc2                                     0.4360736 0.14145323 2554
## Mother.Acc2                                     0.2226406 0.10883690 2554
## bSpline:Offspring.Acc2                         -0.0028834 0.12384042 2554
## bSpline:Father.Acc2                            -0.3593654 0.13311546 2554
## Offspring.Acc2:Father.Acc2                     -0.0476710 0.15979645 2554
## bSpline:Mother.Acc2                            -0.1737977 0.12144948 2554
## Offspring.Acc2:Mother.Acc2                      0.0112290 0.15280483 2554
## Father.Acc2:Mother.Acc2                        -0.2187508 0.15867081 2554
## bSpline:Offspring.Acc2:Father.Acc2             -0.0482553 0.17574242 2554
## bSpline:Offspring.Acc2:Mother.Acc2             -0.1292361 0.16544425 2554
## bSpline:Father.Acc2:Mother.Acc2                 0.1691838 0.18017240 2554
## Offspring.Acc2:Father.Acc2:Mother.Acc2         -0.1522449 0.21895744 2554
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2  0.3038288 0.23939520 2554
##                                                  t-value p-value
## (Intercept)                                    -7.483075  0.0000
## Mass                                           12.272296  0.0000
## bSpline                                        13.026686  0.0000
## Offspring.Acc2                                  0.150528  0.8804
## Father.Acc2                                     3.082811  0.0021
## Mother.Acc2                                     2.045635  0.0409
## bSpline:Offspring.Acc2                         -0.023283  0.9814
## bSpline:Father.Acc2                            -2.699652  0.0070
## Offspring.Acc2:Father.Acc2                     -0.298323  0.7655
## bSpline:Mother.Acc2                            -1.431029  0.1525
## Offspring.Acc2:Mother.Acc2                      0.073486  0.9414
## Father.Acc2:Mother.Acc2                        -1.378646  0.1681
## bSpline:Offspring.Acc2:Father.Acc2             -0.274580  0.7837
## bSpline:Offspring.Acc2:Mother.Acc2             -0.781146  0.4348
## bSpline:Father.Acc2:Mother.Acc2                 0.939010  0.3478
## Offspring.Acc2:Father.Acc2:Mother.Acc2         -0.695317  0.4869
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2  1.269152  0.2045
##  Correlation: 
##                                                (Intr) Mass   bSplin Off.A2
## Mass                                           -0.475                     
## bSpline                                        -0.618 -0.002              
## Offspring.Acc2                                 -0.487 -0.081  0.607       
## Father.Acc2                                    -0.647  0.095  0.494  0.397
## Mother.Acc2                                    -0.559  0.029  0.641  0.532
## bSpline:Offspring.Acc2                          0.460 -0.006 -0.740 -0.841
## bSpline:Father.Acc2                             0.425  0.005 -0.690 -0.420
## Offspring.Acc2:Father.Acc2                      0.365  0.000 -0.438 -0.694
## bSpline:Mother.Acc2                             0.466  0.001 -0.754 -0.458
## Offspring.Acc2:Mother.Acc2                      0.382  0.011 -0.456 -0.735
## Father.Acc2:Mother.Acc2                         0.367 -0.006 -0.441 -0.357
## bSpline:Offspring.Acc2:Father.Acc2             -0.325  0.002  0.522  0.596
## bSpline:Offspring.Acc2:Mother.Acc2             -0.347  0.009  0.554  0.629
## bSpline:Father.Acc2:Mother.Acc2                -0.314 -0.004  0.510  0.310
## Offspring.Acc2:Father.Acc2:Mother.Acc2         -0.252 -0.027  0.321  0.507
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2  0.240 -0.004 -0.385 -0.437
##                                                Fth.A2 Mth.A2 bSp:O.A2 bSp:F.A2
## Mass                                                                          
## bSpline                                                                       
## Offspring.Acc2                                                                
## Father.Acc2                                                                   
## Mother.Acc2                                     0.425                         
## bSpline:Offspring.Acc2                         -0.366 -0.474                  
## bSpline:Father.Acc2                            -0.714 -0.442  0.510           
## Offspring.Acc2:Father.Acc2                     -0.611 -0.373  0.607    0.632  
## bSpline:Mother.Acc2                            -0.372 -0.846  0.558    0.519  
## Offspring.Acc2:Mother.Acc2                     -0.298 -0.710  0.633    0.315  
## Father.Acc2:Mother.Acc2                        -0.608 -0.672  0.326    0.636  
## bSpline:Offspring.Acc2:Father.Acc2              0.541  0.336 -0.707   -0.757  
## bSpline:Offspring.Acc2:Mother.Acc2              0.274  0.621 -0.749   -0.381  
## bSpline:Father.Acc2:Mother.Acc2                 0.527  0.572 -0.377   -0.738  
## Offspring.Acc2:Father.Acc2:Mother.Acc2          0.438  0.486 -0.444   -0.462  
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.397 -0.431  0.519    0.556  
##                                                Of.A2:F.A2 bS:M.A O.A2:M F.A2:M
## Mass                                                                          
## bSpline                                                                       
## Offspring.Acc2                                                                
## Father.Acc2                                                                   
## Mother.Acc2                                                                   
## bSpline:Offspring.Acc2                                                        
## bSpline:Father.Acc2                                                           
## Offspring.Acc2:Father.Acc2                                                    
## bSpline:Mother.Acc2                             0.330                         
## Offspring.Acc2:Mother.Acc2                      0.517      0.602              
## Father.Acc2:Mother.Acc2                         0.537      0.581  0.478       
## bSpline:Offspring.Acc2:Father.Acc2             -0.857     -0.393 -0.448 -0.482
## bSpline:Offspring.Acc2:Mother.Acc2             -0.453     -0.734 -0.844 -0.426
## bSpline:Father.Acc2:Mother.Acc2                -0.467     -0.675 -0.407 -0.859
## Offspring.Acc2:Father.Acc2:Mother.Acc2         -0.727     -0.422 -0.692 -0.725
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2  0.629      0.509  0.586  0.647
##                                                bSp:O.A2:F.A2 bS:O.A2:M bS:F.A2:
## Mass                                                                           
## bSpline                                                                        
## Offspring.Acc2                                                                 
## Father.Acc2                                                                    
## Mother.Acc2                                                                    
## bSpline:Offspring.Acc2                                                         
## bSpline:Father.Acc2                                                            
## Offspring.Acc2:Father.Acc2                                                     
## bSpline:Mother.Acc2                                                            
## Offspring.Acc2:Mother.Acc2                                                     
## Father.Acc2:Mother.Acc2                                                        
## bSpline:Offspring.Acc2:Father.Acc2                                             
## bSpline:Offspring.Acc2:Mother.Acc2              0.528                          
## bSpline:Father.Acc2:Mother.Acc2                 0.559         0.495            
## Offspring.Acc2:Father.Acc2:Mother.Acc2          0.626         0.590     0.623  
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.734        -0.693    -0.753  
##                                                O.A2:F.A2:
## Mass                                                     
## bSpline                                                  
## Offspring.Acc2                                           
## Father.Acc2                                              
## Mother.Acc2                                              
## bSpline:Offspring.Acc2                                   
## bSpline:Father.Acc2                                      
## Offspring.Acc2:Father.Acc2                               
## bSpline:Mother.Acc2                                      
## Offspring.Acc2:Mother.Acc2                               
## Father.Acc2:Mother.Acc2                                  
## bSpline:Offspring.Acc2:Father.Acc2                       
## bSpline:Offspring.Acc2:Mother.Acc2                       
## bSpline:Father.Acc2:Mother.Acc2                          
## Offspring.Acc2:Father.Acc2:Mother.Acc2                   
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.849    
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.02357916 -0.61273768 -0.06675436  0.56340609  3.77352550 
## 
## Number of Observations: 2571
## Number of Groups: 1
# Model output with corrected degrees of freedom and p-values

print(coef_output)
##                                                       Value  Std.Error
## (Intercept)                                    -0.844786431 0.11289296
## Mass                                            0.204952884 0.01670045
## bSpline                                         1.193208038 0.09159720
## Offspring.Acc2                                  0.017313773 0.11502032
## Father.Acc2                                     0.436073634 0.14145323
## Mother.Acc2                                     0.222640621 0.10883690
## bSpline:Offspring.Acc2                         -0.002883363 0.12384042
## bSpline:Father.Acc2                            -0.359365393 0.13311546
## Offspring.Acc2:Father.Acc2                     -0.047671022 0.15979645
## bSpline:Mother.Acc2                            -0.173797699 0.12144948
## Offspring.Acc2:Mother.Acc2                      0.011229015 0.15280483
## Father.Acc2:Mother.Acc2                        -0.218750846 0.15867081
## bSpline:Offspring.Acc2:Father.Acc2             -0.048255289 0.17574242
## bSpline:Offspring.Acc2:Mother.Acc2             -0.129236061 0.16544425
## bSpline:Father.Acc2:Mother.Acc2                 0.169183759 0.18017240
## Offspring.Acc2:Father.Acc2:Mother.Acc2         -0.152244934 0.21895744
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2  0.303828827 0.23939520
##                                                        DF     t.value
## (Intercept)                                      19.04531 -7.48307473
## Mass                                            160.99689 12.27229627
## bSpline                                        2338.65930 13.02668634
## Offspring.Acc2                                  415.63258  0.15052795
## Father.Acc2                                      11.41576  3.08281137
## Mother.Acc2                                     213.08059  2.04563541
## bSpline:Offspring.Acc2                         2340.43281 -0.02328290
## bSpline:Father.Acc2                            2338.81306 -2.69965173
## Offspring.Acc2:Father.Acc2                      767.14285 -0.29832341
## bSpline:Mother.Acc2                            2337.24613 -1.43102876
## Offspring.Acc2:Mother.Acc2                      470.20003  0.07348599
## Father.Acc2:Mother.Acc2                         834.12226 -1.37864582
## bSpline:Offspring.Acc2:Father.Acc2             2343.31538 -0.27457964
## bSpline:Offspring.Acc2:Mother.Acc2             2345.15620 -0.78114569
## bSpline:Father.Acc2:Mother.Acc2                2338.46722  0.93901039
## Offspring.Acc2:Father.Acc2:Mother.Acc2          557.56944 -0.69531748
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 2344.32382  1.26915172
##                                                     p.value
## (Intercept)                                    4.387906e-07
## Mass                                           7.360375e-25
## bSpline                                        1.684939e-37
## Offspring.Acc2                                 8.804211e-01
## Father.Acc2                                    1.000747e-02
## Mother.Acc2                                    4.202078e-02
## bSpline:Offspring.Acc2                         9.814266e-01
## bSpline:Father.Acc2                            6.991185e-03
## Offspring.Acc2:Father.Acc2                     7.655371e-01
## bSpline:Mother.Acc2                            1.525557e-01
## Offspring.Acc2:Mother.Acc2                     9.414506e-01
## Father.Acc2:Mother.Acc2                        1.683736e-01
## bSpline:Offspring.Acc2:Father.Acc2             7.836634e-01
## bSpline:Offspring.Acc2:Mother.Acc2             4.347957e-01
## bSpline:Father.Acc2:Mother.Acc2                3.478224e-01
## Offspring.Acc2:Father.Acc2:Mother.Acc2         4.871460e-01
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 2.045129e-01
summary(ModelTMB.AR)
##  Family: gaussian  ( identity )
## Formula:          
## MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc +  
##     (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 |  
##     Degree.Day) + ar1(Time + 0 | AR_Group)
## Data: BTChallenge_Clean
## 
##      AIC      BIC   logLik deviance df.resid 
##   -993.6   -853.2    520.8  -1041.6     2542 
## 
## Random effects:
## 
## Conditional model:
##  Groups     Name        Variance  Std.Dev.  Corr      
##  Father.ID  (Intercept) 1.977e-03 4.446e-02           
##  Mother.ID  (Intercept) 1.530e-13 3.911e-07           
##  Subject.ID (Intercept) 3.160e-10 1.778e-05           
##  Degree.Day (Intercept) 2.000e-03 4.472e-02           
##  AR_Group   Time1       4.335e-02 2.082e-01 0.89 (ar1)
##  Residual               1.998e-02 1.414e-01           
## Number of obs: 2566, groups:  
## Father.ID, 4; Mother.ID, 4; Subject.ID, 229; Degree.Day, 30; AR_Group, 229
## 
## Dispersion estimate for gaussian family (sigma^2): 0.02 
## 
## Conditional model:
##                                                Estimate Std. Error z value
## (Intercept)                                    -0.78843    0.11501  -6.855
## Mass                                            0.20446    0.01690  12.096
## bSpline                                         1.12275    0.10919  10.283
## Offspring.Acc2                                 -0.03014    0.13249  -0.227
## Father.Acc2                                     0.37947    0.14590   2.601
## Mother.Acc2                                     0.11575    0.12572   0.921
## bSpline:Offspring.Acc2                          0.03452    0.14663   0.235
## bSpline:Father.Acc2                            -0.28727    0.15767  -1.822
## Offspring.Acc2:Father.Acc2                     -0.06369    0.18608  -0.342
## bSpline:Mother.Acc2                            -0.01036    0.14333  -0.072
## Offspring.Acc2:Mother.Acc2                      0.06761    0.17570   0.385
## Father.Acc2:Mother.Acc2                        -0.13267    0.18465  -0.718
## bSpline:Offspring.Acc2:Father.Acc2             -0.01277    0.20730  -0.062
## bSpline:Offspring.Acc2:Mother.Acc2             -0.19074    0.19490  -0.979
## bSpline:Father.Acc2:Mother.Acc2                 0.04054    0.21134   0.192
## Offspring.Acc2:Father.Acc2:Mother.Acc2         -0.14540    0.25177  -0.578
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2  0.30242    0.28061   1.078
##                                                Pr(>|z|)    
## (Intercept)                                    7.13e-12 ***
## Mass                                            < 2e-16 ***
## bSpline                                         < 2e-16 ***
## Offspring.Acc2                                   0.8201    
## Father.Acc2                                      0.0093 ** 
## Mother.Acc2                                      0.3572    
## bSpline:Offspring.Acc2                           0.8139    
## bSpline:Father.Acc2                              0.0685 .  
## Offspring.Acc2:Father.Acc2                       0.7322    
## bSpline:Mother.Acc2                              0.9424    
## Offspring.Acc2:Mother.Acc2                       0.7004    
## Father.Acc2:Mother.Acc2                          0.4725    
## bSpline:Offspring.Acc2:Father.Acc2               0.9509    
## bSpline:Offspring.Acc2:Mother.Acc2               0.3277    
## bSpline:Father.Acc2:Mother.Acc2                  0.8479    
## Offspring.Acc2:Father.Acc2:Mother.Acc2           0.5636    
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2   0.2812    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
View(data.frame(
  "Term" = rownames(coef_output),
  "lme_Beta" = coef_output$Value, "TMB_Beta" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 1]),
  "lme_p" = coef_output$p.value, "TMB_p" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 4])
) %>%
  mutate("Diff" = ifelse(lme_p < 0.05 & TMB_p > 0.05, "*",
    ifelse(TMB_p < 0.05 & lme_p > 0.05, "*", " ")
  )))

# Some subtle differences between models
# Now using Kenward-Roger approximation

require(pbkrtest)

get_kr_df <- function(model_object) {
  L <- diag(rep(1, length(fixef(model_object))))
  L <- as.data.frame(L)
  out <- purrr::map_dbl(L, pbkrtest::get_Lb_ddf, object = model_object)
  names(out) <- names(fixef(model_object))
  out
}

KR_DF <- get_kr_df(DF.Model)

p_vals <- c()
for (i in 1:length(t_vals)) {
  if (t_vals[i] > 0) {
    p_vals[i] <- 2 * (pt(t_vals[i], df = as.numeric(KR_DF)[i], lower = FALSE))
  } else if (t_vals[i] < 0) {
    p_vals[i] <- 2 * (pt(t_vals[i], df = as.numeric(KR_DF)[i], lower = TRUE))
  }
}

DF_coef_output <- data.frame(summary(Model.AR)$tTab)
DF_coef_output$p.value <- p_vals
DF_coef_output$DF <- KR_DF

View(data.frame(
  "Term" = rownames(coef_output),
  "lme_BetaSat" = coef_output$Value, "lme_BetaKR" = DF_coef_output$Value,
  "TMB_Beta" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 1]),
  "lme_pSat" = coef_output$p.value, "lme_pKR" = DF_coef_output$p.value,
  "TMB_p" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 4])
) %>%
  mutate("Diff" = ifelse(lme_pKR < 0.05 & TMB_p > 0.05, "*",
    ifelse(TMB_p < 0.05 & lme_pKR > 0.05, "*", " ")
  )))

# Satherwaitte and KR outputs are very similar - probably owing to such large DF estimates.
# Estimating dfs according to sample sizes of offspring alone.

p_vals <- c()

for (i in 1:length(t_vals)) {
  if (t_vals[i] > 0) {
    p_vals[i] <- 2 * (pt(t_vals[i], df = 228, lower = FALSE))
  } else if (t_vals[i] < 0) {
    p_vals[i] <- 2 * (pt(t_vals[i], df = 228, lower = TRUE))
  }
}

Simple_coef_output <- data.frame(summary(Model.AR)$tTab)
Simple_coef_output$p.value <- p_vals
Simple_coef_output$DF <- 228

View(data.frame(
  "Term" = rownames(coef_output),
  "lme_Beta" = coef_output$Value,
  "TMB_Beta" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 1]),
  "lme_pSat" = coef_output$p.value, "lme_pKR" = DF_coef_output$p.value,
  "lme_pSimple" = Simple_coef_output$p.value,
  "TMB_p" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 4])
) %>%
  mutate("Diff" = ifelse(lme_pSimple < 0.05 & lme_pKR > 0.05, "*",
    ifelse(lme_pKR < 0.05 & lme_pSimple > 0.05, "*", " ")
  )))

# Printing off data frame

{
  T1 <- coef_output[, c("Value", "DF", "t.value", "p.value")] %>%
    mutate("Method" = "Satterthwaite", "Term" = rownames(coef_output)) %>%
    rename(
      "Test_Statistic" = "t.value", "Estimate" = "Value",
      "p" = "p.value"
    ) %>%
    select(Term, Estimate, Method, DF, Test_Statistic, p)

  T2 <- DF_coef_output[, c("Value", "DF", "t.value", "p.value")] %>%
    mutate("Method" = "Kenward-Roger", "Term" = rownames(coef_output)) %>%
    rename(
      "Test_Statistic" = "t.value", "Estimate" = "Value",
      "p" = "p.value"
    ) %>%
    select(Term, Estimate, Method, DF, Test_Statistic, p)

  T3 <- as.data.frame(summary(ModelTMB.AR)$coefficient$cond)[, c("Estimate", "z value", "Pr(>|z|)")] %>%
    mutate(
      "Method" = "z-value", "DF" = NA,
      "Term" = rownames(as.data.frame(summary(ModelTMB.AR)$coefficient$cond))
    ) %>%
    rename(
      "Test_Statistic" = "z value",
      "p" = "Pr(>|z|)"
    ) %>%
    select(Term, Estimate, Method, DF, Test_Statistic, p)

  T4 <- Simple_coef_output[, c("Value", "DF", "t.value", "p.value")] %>%
    mutate("Method" = "n-1", "Term" = rownames(coef_output)) %>%
    rename(
      "Test_Statistic" = "t.value", "Estimate" = "Value",
      "p" = "p.value"
    ) %>%
    select(Term, Estimate, Method, DF, Test_Statistic, p)
}

All <- rbind(T1, T2, T3, T4) %>%
  mutate(Method = factor(Method,
    levels = c("Satterthwaite", "Kenward-Roger", "n-1", "z-value")
  )) %>%
  arrange(Term, Method) %>%
  mutate("p < alpha" = ifelse(p < 0.05, "*", ""))

write.csv(All, "Model_Comparison_Final.csv", row.names = FALSE)

Calculating marginal and conditional r-squared values

r.squaredGLMM(Model.AR)
##            R2m       R2c
## [1,] 0.3932383 0.6577368
# Quite low in comparison to models in Penney et al (In Press), but somewhat reasonable.

Producing marginal mean plots for lme model

require("emmeans")

{
  EM_Means_Off <- emmeans::emmip(Model.AR, ~ bSpline | Offspring.Acc, cov.reduce = FALSE, nesting = NULL, nesting.order = FALSE)
  EM_Means_Mat <- emmeans::emmip(Model.AR, ~ bSpline | Mother.Acc, cov.reduce = FALSE)
  EM_Means_Pat <- emmeans::emmip(Model.AR, ~ bSpline | Father.Acc ~ bSpline, cov.reduce = FALSE)
  EM_Means_Off_Pat <- emmeans::emmip(Model.AR, Offspring.Acc ~ bSpline | Father.Acc, cov.reduce = FALSE)
  EM_Means_Off_Mat <- emmeans::emmip(Model.AR, Offspring.Acc ~ bSpline | Mother.Acc, cov.reduce = FALSE)
  EM_Means_Mass <- emmeans::emmip(Model.AR, ~Mass, at = list(Mass = seq(min(Dat$Mass, na.rm = T), max(Dat$Mass, na.rm = T), 0.1)), CIs = TRUE, type = "response", data = Dat)
}

Raw <- as.data.frame(emmeans::emmeans(Model.AR, specs = "bSpline", by = "Offspring.Acc", cov.reduce = FALSE))
bSpline_vals <- Raw$bSpline
Challenge_Temp <- c()

for (i in 1:length(bSpline_vals)) {
  Challenge_Temp[i] <- BTChallenge[c(which(BTChallenge$bSpline == bSpline_vals[i])), 9]
}

# Removing replicates

Challenge_Temp <- unique(Challenge_Temp)
EM_Off <- as.data.frame(EM_Means_Off$data)
EM_Off$Challenge <- sort(rep(Challenge_Temp, 2))

# Removing prodictions for warm acclimated offspring at cold temperatures

to_rm <- c(which(EM_Off$Offspring.Acc == "2" & EM_Off$Challenge < 19))
EM_Off <- EM_Off[-to_rm, ]
EM_Off$Upper.CI <- EM_Off$yvar + EM_Off$SE * 1.96
EM_Off$Lower.CI <- EM_Off$yvar - EM_Off$SE * 1.96

EM_Plot_Off_Ribbon <- ggplot(EM_Off, aes(x = Challenge, y = yvar, fill = Offspring.Acc, linetype = Offspring.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(fill = "Offspring Acclimation \n Temperature", linetype = "Offspring Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Mat <- as.data.frame(EM_Means_Mat$data)
EM_Mat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Mat$Upper.CI <- EM_Mat$yvar + EM_Mat$SE * 1.96
EM_Mat$Lower.CI <- EM_Mat$yvar - EM_Mat$SE * 1.96

EM_Plot_Mat_Ribbon <- ggplot(EM_Mat, aes(x = Challenge, y = yvar, linetype = Mother.Acc, fill = Mother.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(fill = "Mother Acclimation \n Temperature", linetype = "Mother Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Pat <- as.data.frame(EM_Means_Pat$data)
EM_Pat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Pat$Upper.CI <- EM_Pat$yvar + EM_Pat$SE * 1.96
EM_Pat$Lower.CI <- EM_Pat$yvar - EM_Pat$SE * 1.96

EM_Plot_Pat_Ribbon <- ggplot(EM_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(fill = "Father Acclimation \n Temperature", linetype = "Father Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

Off_Pat <- as.data.frame(EM_Means_Off_Pat$data)
Off_Pat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Pat$Upper.CI <- Off_Pat$yvar + Off_Pat$SE * 1.96
Off_Pat$Lower.CI <- Off_Pat$yvar - Off_Pat$SE * 1.96

to_rm <- c(which(Off_Pat$Offspring.Acc == "2" & Off_Pat$Challenge < 19))
Off_Pat <- Off_Pat[-to_rm, ]

facet_names <- c(`1` = "Cold Fathers", `2` = "Warm Fathers")
facet_names_new <- c(`1` = "Cold Offspring", `2` = "Warm Offspring")

EM_Plot_Off_Pat_Ribbon <- ggplot(Off_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(linetype = "Father Acclimation \n Temperature", fill = "Father Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
  xlim(c(15, 30))

Off_Mat <- as.data.frame(EM_Means_Off_Mat$data)
Off_Mat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Mat$Upper.CI <- Off_Mat$yvar + Off_Mat$SE * 1.96
Off_Mat$Lower.CI <- Off_Mat$yvar - Off_Mat$SE * 1.96

to_rm <- c(which(Off_Mat$Offspring.Acc == "2" & Off_Mat$Challenge < 19))
Off_Mat <- Off_Mat[-to_rm, ]
facet_names <- c(`1` = "Cold Mothers", `2` = "Warm Mothers")

EM_Plot_Off_Mat_Ribbon <- ggplot(Off_Mat, aes(x = Challenge, y = yvar, fill = Mother.Acc, linetype = Mother.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(linetype = "Mother Acclimation \n Temperature", fill = "Mother Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Mass <- as.data.frame(EM_Means_Mass$data)
EM_Mass$Upper.CI <- EM_Mass$yvar + EM_Mass$SE * 1.96
EM_Mass$Lower.CI <- EM_Mass$yvar - EM_Mass$SE * 1.96

EM_Plot_Mass <- ggplot(EM_Mass, aes(x = Mass, y = yvar)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Mass, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  xlab("Mass (g)") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Mass (g)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank())

grid.arrange(EM_Plot_Off_Ribbon, EM_Plot_Mat_Ribbon, EM_Plot_Pat_Ribbon, EM_Plot_Off_Mat_Ribbon, EM_Plot_Off_Pat_Ribbon, EM_Plot_Mass, nrow = 3)

Spline_plot <- BTChallenge %>%
  ggplot(aes(x = Challenge.Temp, y = MO2.NMS)) +
  geom_point(pch = 21, size = 3, alpha = 0.3, colour = "black", fill = "royalblue3") +
  geom_smooth(
    mapping = aes(x = Challenge.Temp, y = bSpline),
    colour = "black", linetype = "dashed"
  ) +
  labs(x = "Challenge Temperature (°C)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
  theme_bw() +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"))

print(Spline_plot)

# Producing plots to save

dev.new()
EM_Plot_Off_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_alone.jpeg", EM_Plot_Off_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mothers_alone.jpeg", EM_Plot_Mat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Fathers_alone.jpeg", EM_Plot_Pat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Off_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_mothers.jpeg", EM_Plot_Off_Mat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Off_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_fathers.jpeg", EM_Plot_Off_Pat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Mass
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mass_alone.jpeg", EM_Plot_Mass, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
Spline_plot
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Spline_plot.jpeg", Spline_plot, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

## With dots

EM_Plot_Off_Ribbon <- ggplot(EM_Off, aes(x = Challenge, y = yvar, fill = Offspring.Acc, linetype = Offspring.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  labs(fill = "Offspring Acclimation \n Temperature", linetype = "Offspring Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30)) +
  stat_summary(geom = "point", fun = "mean", size = 2, colour = "black", pch = 21, alpha = 0.2, data = Dat, aes(x = Challenge.Temp, y = MO2.NMS))

EM_Plot_Mat_Ribbon <- ggplot(EM_Mat, aes(x = Challenge, y = yvar, colour = Mother.Acc, fill = Mother.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  labs(fill = "Mother Acclimation \n Temperature", colour = "Mother Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30)) +
  stat_summary(geom = "point", fun = "mean", size = 2, colour = "black", pch = 21, alpha = 0.2, data = Dat, aes(x = Challenge.Temp, y = MO2.NMS))

EM_Plot_Pat_Ribbon <- ggplot(EM_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, colour = Father.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  labs(fill = "Father Acclimation \n Temperature", colour = "Father Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

facet_names <- c(`1` = "Cold Fathers", `2` = "Warm Fathers")
facet_names_new <- c(`1` = "Cold Offspring", `2` = "Warm Offspring")

EM_Plot_Off_Pat_Ribbon <- ggplot(Off_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, colour = Father.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  labs(colour = "Father Acclimation \n Temperature", fill = "Father Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
  xlim(c(15, 30))

facet_names <- c(`1` = "Cold Mothers", `2` = "Warm Mothers")

EM_Plot_Off_Mat_Ribbon <- ggplot(Off_Mat, aes(x = Challenge, y = yvar, fill = Mother.Acc, colour = Mother.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  labs(colour = "Mother Acclimation \n Temperature", fill = "Mother Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Plot_Mass <- ggplot(EM_Mass, aes(x = Mass, y = yvar)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Mass, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  xlab("Mass (g)") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Mass (g)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank())

Producing marginal mean plots for glmmTMB model

{
  EM_Means_Off <- emmeans::emmip(ModelTMB.AR, Offspring.Acc ~ bSpline, cov.reduce = FALSE)
  EM_Means_Mat <- emmeans::emmip(ModelTMB.AR, Mother.Acc ~ bSpline, cov.reduce = FALSE)
  EM_Means_Pat <- emmeans::emmip(ModelTMB.AR, Father.Acc ~ bSpline, cov.reduce = FALSE)
  EM_Means_Off_Pat <- emmeans::emmip(ModelTMB.AR, Offspring.Acc ~ bSpline | Father.Acc, cov.reduce = FALSE)
  EM_Means_Off_Mat <- emmeans::emmip(ModelTMB.AR, Offspring.Acc ~ bSpline | Mother.Acc, cov.reduce = FALSE)
  EM_Means_Mass <- emmeans::emmip(ModelTMB.AR, ~Mass, at = list(Mass = seq(min(Dat$Mass, na.rm = T), max(Dat$Mass, na.rm = T), 0.1)), CIs = TRUE, type = "response", data = Dat)
}

Raw <- as.data.frame(emmeans::emmeans(ModelTMB.AR, specs = "bSpline", by = "Offspring.Acc", cov.reduce = FALSE))
bSpline_vals <- Raw$bSpline

Challenge_Temp <- c()
for (i in 1:length(bSpline_vals)) {
  Challenge_Temp[i] <- BTChallenge[c(which(BTChallenge$bSpline == bSpline_vals[i])), 9]
}

# Removing replicates again

Challenge_Temp <- unique(Challenge_Temp)
EM_Off <- as.data.frame(EM_Means_Off$data)
EM_Off$Challenge <- sort(rep(Challenge_Temp, 2))

# Removing prodictions for warm acclimated offspring at cold temperatures

to_rm <- c(which(EM_Off$Offspring.Acc == "2" & EM_Off$Challenge < 19))
EM_Off <- EM_Off[-to_rm, ]
EM_Off$Upper.CI <- EM_Off$yvar + EM_Off$SE * 1.96
EM_Off$Lower.CI <- EM_Off$yvar - EM_Off$SE * 1.96

EM_Plot_Off_Ribbon <- ggplot(EM_Off, aes(x = Challenge, y = yvar, fill = Offspring.Acc, linetype = Offspring.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(fill = "Offspring Acclimation \n Temperature", linetype = "Offspring Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Mat <- as.data.frame(EM_Means_Mat$data)
EM_Mat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Mat$Upper.CI <- EM_Mat$yvar + EM_Mat$SE * 1.96
EM_Mat$Lower.CI <- EM_Mat$yvar - EM_Mat$SE * 1.96

EM_Plot_Mat_Ribbon <- ggplot(EM_Mat, aes(x = Challenge, y = yvar, linetype = Mother.Acc, fill = Mother.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(fill = "Mother Acclimation \n Temperature", linetype = "Mother Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Pat <- as.data.frame(EM_Means_Pat$data)
EM_Pat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Pat$Upper.CI <- EM_Pat$yvar + EM_Pat$SE * 1.96
EM_Pat$Lower.CI <- EM_Pat$yvar - EM_Pat$SE * 1.96

EM_Plot_Pat_Ribbon <- ggplot(EM_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(fill = "Father Acclimation \n Temperature", linetype = "Father Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank()) +
  xlim(c(15, 30))

Off_Pat <- as.data.frame(EM_Means_Off_Pat$data)
Off_Pat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Pat$Upper.CI <- Off_Pat$yvar + Off_Pat$SE * 1.96
Off_Pat$Lower.CI <- Off_Pat$yvar - Off_Pat$SE * 1.96

to_rm <- c(which(Off_Pat$Offspring.Acc == "2" & Off_Pat$Challenge < 19))
Off_Pat <- Off_Pat[-to_rm, ]

facet_names <- c(`1` = "Cold Fathers", `2` = "Warm Fathers")
facet_names_new <- c(`1` = "Cold Offspring", `2` = "Warm Offspring")

EM_Plot_Off_Pat_Ribbon <- ggplot(Off_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(linetype = "Father Acclimation \n Temperature", fill = "Father Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
  xlim(c(15, 30))

Off_Mat <- as.data.frame(EM_Means_Off_Mat$data)
Off_Mat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Mat$Upper.CI <- Off_Mat$yvar + Off_Mat$SE * 1.96
Off_Mat$Lower.CI <- Off_Mat$yvar - Off_Mat$SE * 1.96

to_rm <- c(which(Off_Mat$Offspring.Acc == "2" & Off_Mat$Challenge < 19))
Off_Mat <- Off_Mat[-to_rm, ]

facet_names <- c(`1` = "Cold Mothers", `2` = "Warm Mothers")

EM_Plot_Off_Mat_Ribbon <- ggplot(Off_Mat, aes(x = Challenge, y = yvar, fill = Mother.Acc, linetype = Mother.Acc)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
  scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
  labs(linetype = "Mother Acclimation \n Temperature", fill = "Mother Acclimation \n Temperature") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
  labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
  theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
  xlim(c(15, 30))

EM_Mass <- as.data.frame(EM_Means_Mass$data)
EM_Mass$Upper.CI <- EM_Mass$yvar + EM_Mass$SE * 1.96
EM_Mass$Lower.CI <- EM_Mass$yvar - EM_Mass$SE * 1.96

EM_Plot_Mass <- ggplot(EM_Mass, aes(x = Mass, y = yvar)) +
  theme_bw() +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
  geom_ribbon(mapping = aes(x = Mass, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
  xlab("Mass (g)") +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
  labs(x = "Mass (g)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
  theme(panel.grid.major = element_blank())

grid.arrange(EM_Plot_Off_Ribbon, EM_Plot_Mat_Ribbon, EM_Plot_Pat_Ribbon, EM_Plot_Off_Mat_Ribbon, EM_Plot_Off_Pat_Ribbon, EM_Plot_Mass, nrow = 3)

Spline_plot <- BTChallenge %>%
  ggplot(aes(x = Challenge.Temp, y = MO2.NMS)) +
  geom_point(pch = 21, size = 3, alpha = 0.3, colour = "black", fill = "royalblue3") +
  geom_smooth(
    mapping = aes(x = Challenge.Temp, y = bSpline),
    colour = "black", linetype = "dashed"
  ) +
  labs(x = "Challenge Temperature (°C)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
  theme_bw() +
  theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"))

print(Spline_plot)

dev.new()
EM_Plot_Off_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_alone_z.jpeg", EM_Plot_Off_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mothers_alone_z.jpeg", EM_Plot_Mat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Fathers_alone_z.jpeg", EM_Plot_Pat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Off_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_mothers_z.jpeg", EM_Plot_Off_Mat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Off_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_fathers_z.jpeg", EM_Plot_Off_Pat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
EM_Plot_Mass
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mass_alone_z.jpeg", EM_Plot_Mass, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)

dev.new()
Spline_plot
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Spline_z.jpeg", Spline_plot, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)